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The aim of hybrid methods in simulations is to communicate regions with disparate time and 
length scales. Here, a fluid described at the atomistic level within an inner region P is coupled to an 
' outer region C described by continuum fluid dynamics. The matching of both descriptions of matter 

| is made across an overlapping region and, in general, consists of a two-way coupling scheme (C^P 

, and P— »C) which conveys mass, momentum and energy fluxes. The contribution of the hybrid 

scheme hereby presented is two-fold: first it treats unsteady flows and, more importantly, it handles 
energy exchange between both C and P regions. The implementation of the C^P coupling is tested 
here using steady and unsteady flows with different rates of mass, momentum and energy exchange. 
In particular, relaxing flows described by linear hydrodynamics (transversal and longitudinal waves) 
' are most enlightening as they comprise the whole set of hydrodynamic modes. Applying the hybrid 

coupling scheme after the onset of an initial perturbation, the cell-averaged Fourier components of 
the flow variables in the P region (velocity, density, internal energy, temperature and pressure) evolve 
in excellent agreement with the hydrodynamic trends. It is also shown that the scheme preserves 
the correct rate of entropy production. We discuss some general requirements on the coarse-grained 
^ , length and time scales arising from both the characteristic microscopic and hydrodynamic scales. 

^ ; PACS numbers: 02.70 -c, 47.11+j, 47.10 +g, 68.65-k 
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I. INTRODUCTION 



A wide range of systems with important applications are governed by a fine interplay between the atomistic pro- 
cesses occurring within a small region of the system and the slow dynamics occurring within the bulk. A large list of 
examples arise in complex flows near interfaces (polymers or colloids near surfaces, wetting, drop formation, melting, 
0^ ' crystal growth from a fluid phase, moving interfaces of immiscible fluids or membranes, to name only a few). The 
computational expense of realistic-size simulations of these problems via standard molecular dynamics (MD) is pro- 
hibitive, and such kind of studies require new algorithms which can retain the benefit of the atomistic description of 
\ matter where it is really needed, while treating the bulk of the system by much less costly continuum fluid mechanics 
• methods. 

Several hybrid algorithms of this sort have been proposed in the recent literature. In general, to couple the particle 
region P and the continuum region C, such hybrid schemes use an overlapping region comprised of two buffers C^P 
and P— >C, which account for the two-way transfer of information: from C to P and vice versa (see Fig. 1). While the 
P^C transfer essentially consists of a coarse-graining procedure, at C^P one needs to reconstruct the dynamics of a 
large collection of particles with only the limited prescription from the C region as input. Moreover, in performing this 
reconstruction, the number of unphysical artifacts added (as Maxwell daemons) should be minimized far as possible. 
This task is very complicated and represents in fact the heart of any hybrid scheme. 

Hybrid algorithms for fluids are relatively recent. The elegant method introduced by Garcia et al. [l| for rarefied 
gases couples fluxes arising from a direct simulation Monte Carlo (DSMC) scheme to another region described by 
computational fluid dynamics (CFD). The DSMC is set at the finest grid scale of an adaptive mesh refinement 
hierarchy, while a CFD semi-implicit solver is used at the upper level scales. In passing, we note that the scheme 
may, in principle, be implemented using an (MD-Continuum) liquid description, although in this case the C-solver 
must be completely explicit to avoid having to change the particle's energy in the iterations of the implicit scheme. 

In the case of liquids, the state of the art is relatively less developed due to the complications arising from the 
interparticle forces. A pioneering work by O'Connel and Thompson Q coupled momentum by imposing the local 
continuum velocity at C— >P via a crude constraint Lagrangian dynamics. Hadjiconstantinou and Patera |3| introduced 
a reservoir region to impose boundary conditions on the P region (this reservoir being the equivalent of the C— >P 
domain defined here). While in residence in the reservoir, particles were given at each time step, a velocity drawn from 



* R.Dclgado-Buscalioni@ucl.ac.uk 
IP . V . Coveney @ucl .ac.uk 



2 



a Maxwellian distribution with mean and variance consistent with the velocity and temperature of the C-flow. To 
obtain the boundary condition for the C region the authors used a low order polynomial to smooth the field variables 
derived from P at the P^C region. In order to match the boundary conditions for both the P and C regions, 
Hadjiconstantinou and Patera 3] implemented an iterative scheme (based on the Schwarz alternating method) that is 
suitable for steady incompressible flows. Liao and Yip 4] proposed a sophisticated method (called the thermodynamic 
field estimator) to extract continuum fields from the particle data by means of maximum likehood inference. This 
idea may be used to ameliorate the P— >C coupling when the flow presents large gradients, albeit at a rather large 
computational cost. To transfer momentum on the P region Liao and Yip |4| proposed a new Maxwell daemon, called 
reflecting particle method. A drawback is that the pressure gradient is then an outcome of the simulation, rather than 
an input. Finally, Flekkoy et al. |l| used the idea of coupling-through-fluxes and also implemented mass transfer. 
However, energy transfer was still not allowed and only steady flows were considered. The main purpose of the present 
work is to broaden the scope of such hybrid schemes towards a general description allowing mass, momentum and 
energy coupling in unsteady flows. 

A question of central interest is to decide what kind of information needs to be transfered at C^P and P— *C 
There are essentially two possibilities, to transfer either generalised forces (fluxes of conserved quantities) or the local 
values of the averaged- variables. Both kinds of approaches can be found in the published literature. Here, in the 
context of energy transfer, we show that under unsteady flows it is not sufficient to impose the local C-quantities at 
the boundary of P; instead, it is necessary to couple through fluxes. Another possible benefit of flux-coupling was 
pointed out by Flekkoy et al. |5j who stated that this procedure transcends the problem of working with fluids whose 
constitutive relations may be only partially or incompletely known. Although we agree that the flux-based coupling 
is the correct matching procedure, we show nevertheless that if the transport coefficients at C and P are disparate 
enough, the hybrid scheme fails to couple the time evolution of both domains. Hence, in such cases, the evaluation 
of transport coefficients (using standard microscopic techniques, at least for the range of densities and temperatures 
under study) is an unavoidable requirement for the correct behaviour of the hybrid scheme. 

The rest of the paper proceeds as follows. The equations governing the continuum and particle regions and the 
averaging procedures are presented in Sec. [H] The core of the scheme, describing the C^P coupling for momentum, 
energy and mass fluxes, is presented in Sec IIIII General requirements on the coarse-graining length and time scales 
are discussed in Section IIVI The unsteady flows under which the scheme has been tested (decay of longitudinal and 
transversal waves) are presented in Sec. [V] and in Sec. I VII we discuss the results of these tests as far as the main 
hydrodynamic and thermodynamic variables are concerned. Finally Sec. IVIII is devoted to conclusions. 



II. OVERVIEW AND GEOMETRY OF THE HYBRID COUPLING SCHEME 



The domain decomposition of the hybrid scheme is depicted in Fig. 1. Two regions need to be distinguished: 
the particle region P, and the continuum region C. Region P is composed of an ensemble of particles interacting 
through prescribed interparticle potentials and evolving in time through Newtonian dynamics. In order to illustrate 
the coupling procedure a Lennard-Jones (LJ) fluid will be considered. Within P a number N(t) of particles, located at 
r = {n} (the subscript i denoting the i th particle) interacts through the LJ potential tp(r) — 4e _1 [(cr/r) 12 — (cr/r) 6 ] . 
Each particle has a mass m, velocity Vi and energy ej = \mv1 + E^^r^) (r^ = Yj — r^). Their equations of motion, 

r, = v i; (1) 
v, = Vm=^^ll?, (2) 

CLTij Tij 

are solved via standard molecular dynamics (MD) at time steps Atp ~ 10~ 3 r, where r = (mcr 2 /e) 1 / 2 is the charac- 
teristic time of the LJ potential. Throughout the rest of the paper, all quantities will be expressed in reduced units 
of the LJ potential: t(= 0.45 x KT 13 s), cr(= 3.305 x 10~ 12 cm), e, m{= 6.63 x 10" 23 g) and e/k B {= 119. 18K) for time, 
length, energy, mass and temperature , respectively (the numerical values correspond to argon). 

On the other hand, within the C region the relevant variables are the macroscopic local densities associated with 
the conserved quantities, the number density p(R, t), the energy density pe(R, t) and the momentum density j(R, t) 
(related to the local mean velocity u by j = pu) . In what follows the spatial coordinates of the macroscopic fields 
are denoted by capital letters R, while the the microscopic coordinates are designated by lower case letters. The 
conservation laws for the local densities are 



dp 

m 
m 



-V • pu, (3) 

• (ju + n) , (4) 
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doe 

= -V- (peu + II-u + q), (5) 

where the specific energy e = u 2 /2 + 3T/2 + 0, includes the translational energy, the thermal kinetic energy and 
the potential energy The momentum flow contains contributions from convection ju and the pressure tensor 
II = P1 + t. the latter including the local hydrostatic pressure P(R, t) and the viscous stress tensor, which satisfy 
a Newtonian constitutive relation, as shown by previous previous MD descriptions of the L J-fluid 0, U$ j 

T = -i/|vu+(Vuf-^V-uj-^V'U. (6) 

The energy current includes convection peu, dissipation II • u and conduction q, which can be expressed in terms of 
the local temperature gradients and the thermal conductivity k c , through Fourier's law q = -k c VT(R, i). In order 
to close the above equations it is necessary to know the caloric e(p, T) and thermal equations of state P(p,T) and 
the constitutive relations for the transport coefficients (shear and bulk viscosities and thermal conductivity; 77, £ and 
k c respectively) in terms of a set of independent thermodynamic variables, such as p and T. The equations of state 
for a LJ-fluid were extracted from Johnson et al. and the transport coefficients rj, k c and £ from Heyes and 
Borgelt et al. The variables relevant to the C region are the slower ones. Using any standard continuum fluid 
dynamics solver (e.g. based on a finite volume method) the evolution of the C- variables will be traced at time intervals 
Ate >> Atp and evaluated within cells of volume V\ whose size and location is given by the nodes of a certain mesh, 
{R;}, I — {1, ...,M C }. It will be assumed that the size of the C— >P and P^C regions are the same size as those of 
the cells used in the spatial discretisation of the selected continuum solver, say Vi — (AX) 3 . In general both AX 
and Ate may depend on the type of solver used for the C-region, or on the characteristic length of the particular 
phenomena under study, Nevertheless, various intrinsic constraints on AX and Ate will be mentioned in Sec II VI 



1. Averages 

Averages are needed in order to transfer information from the faster time and shorter length-scale particle dynamics 
to the slower and longer coarse-grained description. In order to deal with unsteady, non-equilibrium scenarios, averages 
need to be local on the slower time scale and in the coarse-grained spatial coordinates. For any particle variable, say 
<&i, we define the following averages: 

1 Nl 

*(R„t) = ^ x>, ( ? ) 

($)(R,,tc) = -KJ^l *(R*,t)dt (8) 

where the summation in Eq. (JJJ is made over the Ni particles inside the cell I. 

The averaging procedure is needed to translate the P and C "languages" to and from each domain. This translation 
is done within the overlapping region, where the two descriptions of matter coexist (see Fig. 1). In particular, within 
the P— >C cells the many degrees of freedom arising from the particle dynamics are coarse-grained to provide boundary 
conditions at the "upper" C-level. As long as the number of degrees of freedom is very much larger at P than at C, 
this operation is rather straightforward and is based on the microscopic derivation of continuum fluid dynamics [jj. 
We adopt the approach advocated by Flekkoy et al. [f| , in making the information transferred from P to C to be the 
coarse-grained particle-fluxes of conserved quantities. These are 

pu-n P c = i (s£ c 77ivA -n PC (9) 
V pc \ / 

n • n PC = \ pc ( (vfjfmViVi - ^>:; V ;' r, ; F, ; ) ^ • npc (10) 

q ■ npc - \ pc ^ hfjfmav, - ^E^ry Vi Fy) \ • n PC (11) 

where Npc is the number of particles inside the P— >C cell and npc is the surface vector shown in Fig. 1. 

By contrast, within the C— > P cells, the particle dynamics must be modified to conform to the averaged-dynamics 
prescribed by the continuum description. In other words, one needs to construct a sort of "generalised boundary 
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condition" for the particle dynamics. As pointed out in all previous papers on the subject [lj, y, liJ |5j , this represents 
the most demanding challenge in that one needs to invent a way to reconstruct the microscopic dynamics of a large 
number of particles, based on only a few properties of the local continuum variables. Moreover, to ensure that the 
effect on the inner P-region is minimised, it is crucial to reduce as much as possible the unphysical artifacts, such as 
Maxwell daemons, that are added to the particle dynamics at C^P. The present work is focused on this problem, 
which lies at the heart of any hybrid scheme. 



III. THE C^P COUPLING 



This part of the hybrid scheme can be alternatively stated as the imposition of generalised (mass, momentum and 
energy) boundary conditions on an MD simulation box. To deal with this task we have coupled the particle region to 
a collection of flows (with explicit analytical solution) , which involve the whole set of conserved quantities exchanged 
(mass, momentum and energy). In this sense, in the present work our C-solver is not numerical but rather analytical. 
In particular, we use the initial (non-equilibrium) state imposed at P to calculate the time-dependent analytical 
solution at C. This C-flow is then imposed on the P region during the rest of the simulation, meaning that (appart 
from the initial state) the hybrid coupling used in the tests presented here works in one direction only (from C to P). 



A. Imposing fluxes under unsteady flows 

Following Flekkoy et al. 5], at C— > P we shall communicate fluxes of conserved quantities. These fluxes correspond 
to mass, momentum and energy transfers through the outer interface of the C— > P cell (the W-surface in Fig. 1). 
Flekkoy et al. obtained these fluxes from the values of the continuum variables at the centre of the control cell, 
x = xo, instead of at the exact position of the C— >P interface, x = xw ■ We have found that it is essential to take into 
account this apparently unimportant technicality when dealing with unsteady scenarios. Let us consider a general 
conservation equation, with 

^ + V-J^ = ^, (12) 

where is the flux of <j> and the source term vanishes, S$ = 0, as in Eqs. Integrating over the control cell 

C^P, 

d_ 

di 



4>dV + / J ■ nds = 0. (13) 
v Js 



For illustration we shall restrict analysis to the one-dimensional (ID) situation depicted in Fig. 1. In this case one 
obtains, 



d f 

— J <j)dV - A3 c/>E ■ n w = -A3^ w ■ n w , 



(14) 



where the subscripts E (east) and W (west) denote that the variables are measured at x — xe and x = Xw respectively. 
The surface vectors, and are shown in Fig. 1, and in Eq. (|14|) use has been made of nw = — ng. The R.H.S. of 
Eq. (|14fl is the flux current of <fi through the interface W of the control cell, which is precisely the (generalised) force 
we want to introduce on the particles at the C^P buffer. We note that only under steady flows does 3<f> w = J^q 
(to see this, integrate Eq. I|12|) from x = xo to x — xw)- hence only in this case does the evaluation of the fluxes at 
xw using the continuum variables at xo lead to the same converged steady state as if the variables at xw were used 
(although the transients may of course differ). It is possible to provide an estimate of the global error arising from 
evaluating the flux at a position xo shifted S AX with respect xw, over a certain time interval At. In the case of the 
momentum equation the deviation of the stress contribution to the momentum flux J = J • n at any instant would be 
of order AJ ~ V J SAX, with S = \xo — xw\/AX is the distance to the C^P interface (5 = 0.5 in Fig. 1). Assuming 
that the mean velocity field can be expressed as u — u • n ~ vS k ' exp(i/ca;) (k being the dominant wavenumber) and 
using the Newtonian constitutive relation for the viscous tensor in Eq. ©, one obtains A J ~ rj^k 2 u^SAX , where 
rjL = 4r;/3 + £ is the longitudinal viscosity. As a particular example we consider a longitudinal wave and evaluate 
the error along a cycle, At = 2ir/(kc s ) (where c s is the sound velocity). As a crude estimate, the accumulated error 
of the cell- aver aged momentum j • n — p e u is of order p e Au ~ A J At/ AX; using p e = 0.5, r/L — 1 and c s ~ 5, one 
obtains Au/u^ ~ 2it5riLk/{c s p e ) ~ 0.3, for the typical wavenumbers considered here k ~ 0.2. Simulations carried 
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out with the momentum flux evaluated at xo yield relative errors of the averaged velocity at C^P of the same order of 
magnitude as this estimate (see Sec I VI C"2"l and Fig. 6b). We observe that most computational fluid dynamics (CFD) 
codes provide the continuum variables at the centre of the control cells (i.e., at xo), so that in order to evaluate the 
fluxes pertaining to the C— >P exchange it would first be necessary to make use of simple interpolation techniques. 

The fluxes arising from the continuum equations (on the R.H.S. of Eqs. (|15I17() ) are imposed on the particle 
ensemble at the C— >P cells through expressions involving atomistic variables (those on the R.H.S. of I|15I17|I ). 

Apu ■ n, (15) 
A (puu + n) • n, (16) 

A(pue + II • u + q) • n. (17) 

where henceforth, n indicates the vector of the outermost interface of the C— >P cell, pointing towards C. The nomen- 
clature used here follows that of Flekkoy et al. jB|: s(t) indicates the number of particles inserted (s > 0) or removed 
(s < 0) from C— >P per unit of time; the velocity of the inserted/removed particles is v'; F ext is the external force 
applied to each particle i within the C^P cell. The total external force is '£ NcF F? xt , where the summation is over 
the Ncp(t) particles inside C^P. Finally, (e') indicates the energy of the inserted/removed particles and (Jq^) refers 
to an externally imposed heat current. 

As mentioned by Flekkoy et al. insertion of Eq. (|15l) into Ea. i|16|) shows that the rates of change of mo- 
mentum due to convection and local stresses are correctly introduced if (v') = u and ^X)f CP F^ xt ^ — = All ■ n = 

—A (Fn + T • n), respectively. Note that Pn is the hydrostatic pressure force (pointing inwards the C— >P cell), 
while the viscous contribution — r • n depends on the local velocity gradient. 

The balance of the energy flux requires some extra conditions. In Eq. I|17fl the convection, dissipation and conduction 

of energy are balanced if (e') = e, (V^f Ff* ' = ATI • u • n and (Jq * ■ n) = iq-n respectively . 

Let us now consider in more detail how the scheme deals with momentum, energy and mass transfer from C to P. 



I Net 



ms 



<v') + ( £ Ft 



INcf 



ms (e 



'> + ( £ Ff* • v, ) - <J§*) • n = - 



B. Momentum exchange 

The condition (v'} = u ensures the balance of momentum convection. If the mass flux points towards the P region 
(s > 0), this condition is fulfilled by choosing the velocity of the inserted particles from a Maxwellian distribution 
according to the local temperature at the C^P cell, P (v') = (l/27rmfcT) 3 / 2 exp (— m(v — u) 2 /2mfcT). Concerning 
particle removal (s < 0), we note that if the average velocity at the C^P cell is equal to the continuum velocity 
(v) = u, then the average velocity of the subset of extracted particles would be precisely (v') = (v) = u. Hence the 
condition of velocity continuity at C^P is needed to ensure the correct balance of momentum convection. 

The change of momentum due to local stresses establishes the overall external force exerted on the particle region. 
Therefore one needs to determine how the overall external force is distributed on each individual particle. Flekkoy 
et al. distributed the force according to a certain function g(x) satisfying g(x\y) = oo, g(xo) = g'(xo) = 0. 
Normalisation leads to 

pext = _ ( gfo) ] Au . (18) 

where x runs perpendicular to the C— >P interface and the applied force is made constant along each Ate- As g(x) 
tends to infinity as x — > xyy, the applied force diverges as one approaches the C^P interface; hence the density 
nearby x = xw is very small or zero. The function g(x) is thus endowed with a two-fold purpose: it ensures a limiting 
extension to P (as the hydrostatic pressure force always points towards the P region, particles will never cross the 
C-^P interface outwards) while also guaranteeing the existence of a small region where particles can be inserted with 
very low risk of overlapping. 

Despite the benefits of the g(x) function for distributing the externally imposed momentum, we decided to use 
g(x) — 1 for all x inside the C^P cell. The reasons for this choice will become clear when explaining the energy 
exchange, below. The first implication of g(x) — 1 is that the external force is equally distributed among all the 
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particles within the C^P cell. In other words, Ff xt no longer depends on the particle label 

F ext = F ext = _ t _J_ \ Au . n ^ 

\NcpJ 

The second implication of g(x) = 1 is that the particle density profile near the C^P interface no longer vanishes, 
so one needs an efficient way to resolve the problem of overlap on the insertion of new particles. This task is carried 
out by the usher algorithm, as explained below. Finally in order to ensure a finite extent of the particle region, if a 
particle (?) crosses outwards the C— >P interface (in Fig. 1, as» = xw — 8, with 5 > 0) with velocity v, it is substituted 
by another one (j) with y.j — yf, Zj = Zi, Xj — xw +^ and with Vj = Vj. In this way, the overall momentum is strictly 
conserved before and after the particle exchange. 



C. Energy exchange 

1. Advection 



The balance of advected energy requires that (e') = e = u 2 /2 + 3T/2 + <j>. Decomposing the particle energy 
into the kinetic and potential parts e' — [v'] 2 /2 + ip', one sees that, since the new particles are drawn from a 
Maxwell distribution, ([i/] 2 /2) = u 2 /2 + 3T/2. By contrast, the balance of potential energy (ip 1 ) — (f> is much less 
straightforward to implement. This condition is fulfilled by the USHER algorithm, described below. 



2. Dissipation 



One needs to satisfy the following balance of heat dissipation: 

= All u n. (20) 




This condition does not generally hold if the external force is distributed according to an arbitrary g(x). Indeed it 
is not even clear that a function g(x) exists satisfying g{xw) — ► oo and enabling the heat dissipation balance in Eq. 
(|20|l . In any case, such a function g would depend on the particle's velocity distribution and then the problem of 
finding g would become a formidable task at each time step. 

The advantage of using g(x) = 1 now becomes clear. As long as Ff* does not depend on the particle label, one 
can greatly simplify the L.H.S. of Eq. l|2"U|l equation to obtain, 




= N C pF ext ■ (v) = -A II • u • n. (21) 

The last equality follows from construction of the overall force NcpF ext = — Att n and from the continuity of velocity 
<v) = u. 



3. Conduction 



The condition (j^*) • n = Aq • n requires the establishment of a heat current along the C— >P cell representing 
the conduction of energy. This may be implemented by various means; for instance, following the idea of Evans 
and coworkers @ one may include an extra force which pulls the "hotter" particles towards the direction of the 
heat flux and conserves the overall momentum. Alternatively, one may try to impose a Chapman-Enskog velocity 
distribution with the desired heat flux, at some region inside the C^P buffer. In this work we have made use of the 
phcnomenological Fourier's law, q = -k c VT. A temperature gradient is imposed along each C— >P cell by using a 
set of Nose- Hoover thermostats (NHT) placed along the direction of the heat flux. The outer and inner thermostats 
are located a distance d apart, and the temperature difference between both set to d VT • n. Typically, at each C^P 
cell we have used a set of two or three NHT's along a distance of 3a or 4er. The values of the Q-parameter appearing 
in the NHT formulation |l0j were chosen small enough to minimise unphysical dynamics, i.e., we have chosen Q ~ 5. 
The main benefit of using the Nose-Hoover formulation is the small distortion these thermostats introduce to the 
particle dynamics compared with other ways of implementing thermostatting |10| . 
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D. Mass exchange: particle insertion 

One important condition on the particle insertion, inherited from the balance of potential energy, is (ip') — <f> (see 
Sec. IIII C|l . To deal with this task, our strategy has been to place the new particles in positions where -0' — <t>- As 
4> is, roughly speaking, the energy needed to insert a new particle, we shall insert particles with energies close to 
the local chemical potential. In any case, this implies that the insertions need to be made in very precise positions 
which depend on the configuration of the rest of the particles. To this end we have developed an algorithm (usher) 
that guides each new inserted particle to a position where the potential energy is equal to <f> (up to a pre-specified 
threshold). A brief explanation is given below (see Delgado-Buscalioni and Coveney for further details). We first 
note that while the USHER algorithm guides a new particle to a correct location, the rest of the particles remain frozen 
in position, usher essentially performs the following steps: 

1. Place the new particle (i = N + 1) at an initial position inside C— >P, r' ^. 

2. Evaluate fjv+i = YljLi fzv+ij an d 




Typically one can use 5r ~ a. 

3. Move the new particle according to the update rule 

r («+i) - -00 + if(«) 
1 — 1 ^ 2 N + 1 

where 8t — min(A t , 8t a ), with A t ~ 0.05 in reduced units 

4. Evaluate the relative difference between the specific internal energy of the new particle, V'jv+i' an d that pre- 
scribed by the continuum, <fi: Err — IV'at+i — ^l/l'/'l 

5. The particle is correctly inserted if Err is small enough (typically ~ 0.1). 

Let us show how the usher algorithm easily overcomes the problem of possible initial overlap with pre-existing 
particles. An overlap leads to very large values of the interparticle force, Jn+i » 1, so in this case 5t a « 1 and 
5t = 5t a . But by construction, during the interval 5t a the new particle moves a distance of the order of the particle size 
(j in the direction of minimum energy, just enough to avoid any initial overlap. Then, as the particle steadily moves 
towards a local minimum of energy, /jv+i decreases, and 5t a increases until it becomes larger than A t . Then, 5t = A t 
is fixed. For liquid densities varying between p = 0.5 — 0.8, the USHER algorithm typically needs 15-90 iterations 
(single-force evaluations) to correctly place a new particle [13] ■ By introducing the particles with xj)' = <p(l ± Err), 
we found that, upon averaging over Ate, the condition (tp'} = <f> holds within about 2% (even using values of Err as 
large as 0.5). 

Until now we have not mentioned any limitation on the sizes of the coarse mesh and time step AX and 5tc- 
Comments on this topic are very scarce in the previous literature on hybrid methods for fluids. Moreover, since the 
local averages are made using these spatial and temporal windows, it is also appropriate to formulate any condition 
on AX and Ate before presenting the results of the tests carried out for different flows. 

IV. LENGTH AND TIME SCALE PREREQUISITES 

1. Arising from consideration of the microscopic description 

Continuum fluid dynamics rest on the local equilibrium assumption. This means that, in order to define a local 
thermodynamic and hydrodynamic state characterising the coarse-grained variables at the C— >P and P— >C cells, the 
size of these cells needs to be greater than the mean free path A extracted from the particle dynamics. Moreover the 
local equilibrium within each cell should be attained on times scales smaller than Ate- This means that Ate has to 
be larger than the collision time r co ; . In summary, AX > A and Ate > T coi ■ In the case of a L J-fluid it is possible to 
use the hard-sphere approximation to make an order-of-magnitude estimate, AX > 0.2p _1 and At > 0.14p" 1 T~ 1 / 2 . 
These conditions become less restrictive at larger densities; as an example, for T = 1, p = 0.5 and a typical MD 
time-step Atp ~ 10~ 3 , local equilibrium would require Ate/Atp > 100 integration steps. 
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2. Arising from consideration of the hydrodynamic description 

Conditions on AX and Ate are firstly imposed by the smallest characteristic length and time scales involved in 
the process under investigation (say 27r/fc max and 27r/w max respectively). Practically, to correctly recover the smallest 
spatio-temporal flow pattern one needs at least 8 points per period, so Ate < ir/4w max and AX < 7r/4fc max . The 
numerical stability of the C-solver algorithm may also impose limitations. As mentioned in Sec. I, algorithms with 
explicit time discretisation are better suited for the C-solver of a hybrid scheme. A necessary condition for their 
numerical stability is C — U Ate / AX < 1/2, where C is the Courant number and U the maximum characteristic 
flow velocity. The value of U depends on the physical process one is dealing with but, to provide numbers in 
the present discussion, let us assume that we are dealing with low or moderate Reynolds numbers. Then, if the 
process is a diffusive one, U = v/AX; alternatively, if sound waves are relevant within the flow, U = c s . In 
summary, the computational window for Ate should be, 0.14p~ 1 T~ 1 / 2 < Ate < AX/(2U). Using the maximum 
grid-spacing allowed, AX — 7r/4fc max , one obtains the computational windows for Ate shown in Fig. 2 versus p, 
for fc max = {0.1,0.2}, U = AX,c s } and T = 2.5. As expected, a sound wave requires smaller time steps than a 
diffusive process. For large enough fc max , the temporal and spatial computational window may be highly localised; 
e.g. for p — 0.5 one should use 0.2 < At c < 0.5 if waves with wavelength smaller than 30ct need to be captured 
by the coarse-grained description. As the density decreases these conditions become much more restrictive, until the 
acoustic time finally becomes smaller than the collision time (see Fig. 2). Also, in rarefied gases, if p < k ma x/^, the 
mean free path becomes larger than the wavelength, but here we shall not be concerned with situations where the 
Navier-Stokes equations are not appropriate (see Garcia et al. for further discussion pj). 



V. TESTS: HYDRODYNAMIC MODES 



As already mentioned our hybrid scheme has been tested under stationary and unsteady flows. Typical stationary 
nonequilibrium states were considered, such as heat conduction profiles 0] and Couette profiles [a, El- The mi- 
croscopic reconstruction of these flows has been well studied in the literature that nothing new may added here. In 
passing, we note the transient times to achieve the steady state from the rest solution were found to be in agreement 
with the diffusive times L x /k and L x /v. The rest of the discussion will be focused on our choice of unsteady scenarios, 
which are described by the decay of transversal and longitudinal waves. These flows are now briefly presented, using 
standard hydrodynamics. 

Consider a fluid at equilibrium characterised by homogeneous mass density, p e , specific energy e e and a vanishing 
mean velocity u e = 0. Our procedure is to perturb this equilibrium state with different hydrodynamic fields {p p , u, e p } 
(periodic in the x direction, i.e. k = fci) and then make use of the C— >P coupling scheme described in Sec. Illll to verify 
that, within the particle region, the subsequent evolution towards equilibrium is carried out in a hydrodynamically 
consistent way. If use is to be made of the linear hydrodynamic theory, the externally induced perturbations should 
be small enough to guarantee that the relaxation process is always governed by the linearised mass, momentum and 
energy equations ©-10 l^j- We defer further discussion of this point to Sec. IVII below. 

As is customary, to solve the linearised set of equations a Laplace- Fourier transform (LFT) is first performed ^4| . 
The LFT of any perturbative variable, say <I>(r,i), will be denoted as 



/oo 
$(r,i) exp(-ik • r)dr, 
-oo 



dz exp(zzi)<E>(k, t)dt. 



(22) 
(23) 



The LFT of the linearised equation 10- © leads to the following algebraic system for 3? = (jP,T p , jx, jy,jz^j HH , 

M$ T (k,z) = $ T (k,0), (24) 



where the hydrodynamic matrix is, 

/ 

M = 



— iz 





ik 





— iz + Kjk 2 


ik^- 




ikD 


—iz + bk 2 





















vk 2 








-iz + vk 2 



\ 



(25) 
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We note that instead of using e p , the energy equation is expressed in terms of temperature fluctuations. Also, for 
clarity it is better to write the solution <& in terms of the t = Fourier-transformed perturbative heat density Q p and 
pressure P p . These quantities are related to pP and T p through the relations 

eP = c * TP+ (!) T /' (26) 

Q p = p e c v (t p - ^-V) , (27) 



p^a 



mc 2 



P p = 1 (p e a T p - (f) . (28) 

7 

Here c p and c v are the specific heat at constant pressure and volume, respectively; 7 = c p /c v , is the adiabatic 
coefficient; n = K c /c p p e is the thermal diffusivity; a = —(dp/dT)p/p e the thermal expansion coefficient; and the 
kinematic longitudinal viscosity, b = (4?//3 + £)/p e is related to the sound attenuation coefficient, T, through 2T — 
b + (7 — 1)k. Finally D — (dP/dT) p and the adiabatic speed of sound is c 2 — 7 (dPdp) T . 

Provided that, in the hydrodynamic limit, the wavelengths of the perturbations are much larger than the mean 
interparticle distance, it is sufficient to obtain the solution of Eq. 124(1 up to 0(k 2 ) 0|. Putting z->w£R into Eqs. 
(|24|l and l|25|) . leads, after some algebra, to the following identities 

f = E K Q(k,0) + J_± EsrP{kq) + ZZ± E Sfix(k , ), (29) 

p e c p p e a p ac s 

p(k,cj) = --E K Q(k,0) + -^E SR P(k,0) + -E S ij x Qi,0), (30) 

j x (k,w) = — S SJ P(k,0) + J B SflJa (k,0) ) (31) 
mc s 

jyfcto) = E v j y (k,0), (32) 

j,(k,w) - ^(k,0), (33) 

and. using Eqs. J23 and 

Q(k,w) = S«Q(k,0), (34) 

P(k,w) = ^ S fl J P(k,0)+mc s £;s/ix(k,0), (35) 

where the following propagators have been introduced, 

E K (k,(j) = exp(-Kfc 2 t), (36) 

Esn(k,ui) = exp (-rfc 2 t)cos(c s M)> (37) 

E SI {k,uj) = -iexp(-rfc 2 t)sin(c s fci), (38) 

E u (k,u) = exp(-^fc 2 f). (39) 

The shear or transverse modes correspond to momentum perturbations along either y and/or z axis (i.e. perpen- 
dicular to the wavevector fci); from Eqs. (|32(l and (|33|l it is clear that they are completely decoupled. The remaining 
hydrodynamic variables {p p ,T p ,j p } are coupled and conform to three longitudinal modes which can be divided into 
two subgroups: two sound modes (involving longitudinal momentum and pressure at constant entropy) and one heat 
mode (involving heat diffusion). The fact that the heat density is an independent mode is evident from Eq. H34fl . 
where it is seen that it relaxes diffusively oc exp(— nk 2 t). We note that Q p = T e s p is essentially the fluctuation of the 
entropy density s p [Tij|. From Eqs. 1)3 111 and (|35|l it is easily shown that the two sound modes (P/c s ±j x ) decay like 
exp(ikc s t - Tk 2 t). 



VI. RESULTS 



A. Set-up and initial states 



The coupling scheme was implemented and tested in the set-up shown in Fig. 3. The system is periodic along y and 
z directions and the gradients of the continuum variables are set along the x direction. The particle region occupies 
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a region of size L x around x = and of size L y = L z along the periodic directions. The P region is divided into 
control cells of size AX where in local averages are taken. The centres of the two C— >P slabs (the outermost cells) 
are situated at x = ±\L X — AX/2\. The deviation from the local equilibrium assumption was monitored in terms of 
the relative difference of the cell-averaged pressure and energy with respect the values given by the equation of state 
of Johnson et al. Ref |6j. Around a distance 1.5a away from the C— >P interface the typical maximum deviations were 
only about 6%. 

The initial perturbative flow was prepared by first letting the P region relax until a vanishing and homogeneous 
mean flow was obtained. Then, during a small time interval (~ 3t), the particle velocities were periodically changed 
according to a Maxwellian distribution with the desired velocity profile and local cell temperature. The resulting 
initial state was then analysed to extract the Fourier components of the whole set of flow variables (v, p, T, e, P). 
For the sake of consistency these were extracted by Fourier transforms of the cell-averaged variables 

M c 

^£ 4>{X h t) cos(fc„AU (40) 
c i 

&) = ^E^'*) sin (")' ( 41 ) 
c i 

where k n — nk (n £ AT) and; c n = 1 for n — 0, c„ = 2 otherwise. In any case, it was checked that the Fourier 
transform of the microscopic variables cj>^ = c n J2^ <fi(xi,t) exp(—ik n Xi)/N yields essentially the same output as 
Eqs. gDJ} and gTJ. 

The initial Fourier transforms calculated from Eqs. (|40|l were injected into Eqs. I|29l) - (|35[) to obtain the time- 
evolution of the continuum variables. These, in turn, were used to calculate the fluxes imposed on the C^P cells over 
time. The transport coefficients used were those reported in the literature 0,0- As an interesting check, it was found 
that the coupling scheme failed significantly if the transport coefficients used in the C-rcgion differed by more than 
about 15% from those of the LJ fluid. In particular, the oscillation frequency of the averaged particle velocity differed 
with respect to that of C, while correlations decayed at a faster rate than those of the C-flow. Therefore, in cases 
where the constitutive relations are not known, this result means that it is first necessary to measure the transport 
coefficients from the particle dynamics (using any standard molecular technique), before applying the hybrid scheme, 
particularly if unsteady flows are to be studied. 

The wavelengths of the initial perturbations were chosen to be much larger than the mean free path, i.e. 2n\/k << 1, 
in order to work within the hydrodynamic regime. In other words, the dependence of the transport coefficients on the 
wavenumber was negligible |l4| . The amplitudes of the initial perturbation were chosen small enough to ensure that 
the subsequent relaxation process could be described by linear theory. In particular, if v^ 1 ) (t) is the maximum Fourier 
amplitude of the velocity, the typical values of the Reynolds number at t = [Rc= jv^ \p e /(kr])] were Re(0) < 3. As 
| v^ 1 ^ (t) | decays exponentially, convection was present only in the first stages of the relaxing flow, but it was not strong 
enough to produce significant deviations from the linear theory (non-linear effects become dominant for Re> O(10) 
|l5J). The maximum Mach number was less than 0.2, and density fluctuations were around \p (1) \/p e ~ Ma 2 ~ 0.05. 

In another test the hybrid scheme was applied to a fluid in mechanical (u = 0) and thermodynamical equilibrium 
(p = 0.5, T = 3.5, e = 2.7, P — 3.2) during a longer simulation (50r) to check for any possible spurious drift in 
the overall momentum and energy (note that Eq. (|15fl ensures the mass conservation by construction). During this 
calculation, the total momentum inside the P region was conserved up to 5 x 10~ 4 and the total energy fluctuated 
~ 5% around its equilibrium value. The size of these fluctuations is consistent with the system size (which contained 
iV = 1600 particles and a specific heat of c v = 1.8). Note that the total energy of the system cannot be conserved 
because a part the system is connected to a thermostat and it also receives mechanical energy from C. 

B. Transversal waves 

In order to test the transfer of momentum flux along the direction perpendicular to the C^P interface, planar 
shear waves along x were excited in a LJ-fluid with p — 0.5, T — 2.5 and r\ = 0.75 ± 0.05 (the error bar comes from 
Ref.Q). These waves were created by imposing sinusoidal y- velocity profiles v y (x) = Vy 1 ^ sm(kx). 

Figure 4 shows results for perturbations with Vy^ in = 1.0 in a system with L x — 20a, L y = L z = 7a and M c = 10 

cells (AX = 2). The filled circles in Figure 4a correspond to the main Fourier component of the velocity v^l in for 
a calculation using wavenumber k — 0.31. The relaxation process is indeed exponential, until the noise-to-signal 
ratio becomes large enough (around t > 10) and, the observed decay rate (0.14 ± 0.1)r _1 agrees perfectly with the 
theoretical value k 2 rj/p e . To explicitly appreciate the effect of the hybrid coupling we show (open circles in Fig. 4a) 
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the outcome of a calculation in which the unsteady shear stress contribution to the external force F ext was set to zero, 
leaving just the contribution of the equilibrium hydrostatic pressure. As expected, if less momentum flux is provided 
towards P, the decay rate is appreciably slower (0.07r _1 ) than the hydrodynamic one (0.14r _1 ). Figure 4b shows the 
autocorrelation function (ACF) of v9^ n (t) for another perturbation with an slightly different wavenumber, k = 0.35. 
The good agreement with the theoretical decay [exp(— 0.17f), in dashed line] shows that the scheme is able to deal 
with small variations in the perturbation shape. 



C. Longitudinal perturbations 

Longitudinal waves transport mass, momentum and both mechanical and thermal energy; they are therefore per- 
fectly suited for an overall test of the hybrid scheme behaviour. In our set-up the wavevector of these waves was 
perpendicular to the C^P interface and they were generated by imposing mean velocities along the x direction as 
either pure cosinusoidal, sinusoidal profiles, or combinations of both profiles. As explained before, the particle veloci- 
ties were extracted from Maxwell-Boltzmann distributions at the local mean velocities and at a constant temperature. 
The mean velocity profile induces pressure, density and energy fluctuations with a periodic pattern in the x direction. 
The temporal behaviour of the main flow variables is now described. 



1. Mass 



One of the main problems involved in the hybrid mass transfer at C— *P is that although the continuum flux 
is a floating point number, one can only possibly exchange an integer number of particles. In order to adhere 
closely as possible to the prescribed continuum mass flux, the following procedure is followed during each interval 
tc < t < tc + Ate 1 , tc — mAtc, m 6 M. First one evaluates the quantity 



tc+Atc 



s(t)dt. 



(42) 



tc 



This floating point number, which represents the number of particles that should cross the C— >P interface along 
tc < t < tc + Ate, is converted into an integer 5N(t) by the following construction: 



5N(t) = NINT[£(i c ) 

S£(t - t k ) = INT 



m - t k ) 

t-Atr 



(Z(t')-5N(t'))dt' 



(43) 
(44) 



where t k is such that \S£(t — tk)\ < 1 and = t n < tk < ifc+i. The deviation 6£(t — t k ) assimilates the errors made 
through successive rounding off (£ — > NINT[£]). When \5£(t — t k )\ becomes larger than one (at t = {t k }) a particle is 
added to (or extracted from) 6N and the corrector 6£ is then reset to zero. To minimise the effect on the remaining 
particles over each interval Ate, the particle crossings are regularly separated in time, at a rate as close as possible 
to 5N(tc)/ Ate- As illustrated in Fig. 5, this kind of procedure enables us to follow rather closely the desired mass 
flux. 



2. Momentum 



Figures 7, 8 and 9 show the time dependence of the Fourier components of the main hydrodynamic and thermody- 
namic variables. Dashed lines correspond to the theoretical trends obtained via the inverse Fourier transform of Eqs. 

For the reasons explained in Sec lIII A1 it is particularly important to ascertain that the averaged velocity at C^P 
correctly follows the desired continuum flow velocity. Figure 6 shows the instantaneous v x and time averaged velocity 
(v x ) at both C— >P cells. We obtained very good agreement, the deviation of (v x ) with respect the continuum value 
being less than about 5% along most part of the damped oscillation (see Fig. 6). The importance of evaluating the 
continuum fluxes at precisely the C— >P interface (xw in Fig. 1), is illustrated in Fig. 6 by comparison with the 
outcome of a calculation in which the fluxes are evaluated at x = xq. For this calculation, the deviation with respect 
to the continuum prescription is clear and its magnitude agrees with the estimate made in Sec. IIII Al ~ 30%. 

Figure 7a shows the time evolution of the Fourier amplitudes of the velocity in one of the calculations (a pure 
cosinusoidal perturbation with i>i*cos = 0.6). The overall velocity wi°cos remains close to zero, confirming that the 



12 



method conserves the initial total momentum. The v x '■ component (initially set to zero) has been also included 

to display the level of noise. The autocorrelation function (ACF) of Vx — V( x cos ) + W(^ S in) ^ s snown m Fig- 7b 
for two runs with different initial profiles (for more details see the caption of Fig. 7). These data were fitted to 
the hydrodynamic expression arising from Eq. 1)31(1 and in Fig. 7 these fits are shown with dashed lines. The best 
fit to the velocity was obtained with Tk 2 = 0.076 and c s k = 0.867, while for the ACF it yielded Tk 2 = 0.072 and 
c s k = 0.867. These values coincide, within the error bars, with those imposed by the C-flow (obtained upon insertion 
of the transport coefficients reported in the literature pj) 0.071 and 0.88, respectively. 

3. Thermodynamic variables 

We start by observing that, unlike longitudinal momentum and pressure, the relaxation of density, temperature and 
energy perturbations includes not only an acoustic part, but also an entropic contribution proportional to exp(— nk 2 t) 
[see Eqs. I|29[) - (I39[) ]. However, as the initial perturbation considered was a mechanical one 

(Tin) _ the entropic 

contribution is rather small. This can be seen in the theoretical expressions written in Fig. 6: the amplitude at t = 
of the entropic part of any heat-related variable is nearly six times smaller than its mechanical counterpart. This 
observation led us to a more careful study of the effect of heat conduction. As explained in Sec IIII C"3l heat currents 
through each C— >P cell were created by imposing several (typically two) Nose-Hoover thermostats (NHT's) placed 
close to the C^P interface over a distance of 4cr. We found it very informative to study the effect of the number p 
of NHT's per C— >P cell (notation: p-NHTcp) on the collective behaviour of the system. In some calculations, this 
number was reduced to merely 1-NHTcp, in such a way that only the local temperature prescribed by the continuum 
was imposed (and not the heat flux). The results of this comparison may be seen in Figs. 8 and 9. Calculations using 

1- NHTcp yielded essentially the same evolution of the velocity and pressure as those with a larger number of NHT's, 
although the pressure for the 1-NHTcp case showed a slight phase-lag, see Fig. 8d. This is not surprising as v and 
P are governed by acoustic terms and are independent of heat transfers. Using 1-NHTcp, deviations on p p , e p and 
T p with respect to the hydrodynamic trend are indeed appreciable, while results with 2-NHTcp adhere closely to the 
analytical curves (see Fig. 8). In any case, the calculation of the entropy production is the best way to highlight the 
completely different qualitative behaviour of 1-NHTcp with respect to two (or more) NHTcp. 

^ . Entropy 

The entropy fluctuations, or more precisely the perturbation of heat Q p = T e s p , was calculated using Eq. I|27|l . 
We note that, in order to reproduce the steady diffusive heat decay in Eq. (|30[1 . an exact cancellation of the acoustic 
oscillations coming from p p and T p is necessary [see Eqs. (I30f) and (1291) ] . As stated above, the contribution associated 
with heat diffusion in Eqs. (|30[) and 129fl is rather small compared to the mechanical one. This observation indicates 
that small mechanically-driven fluctuations around the local thermodynamic equilibrium may become large enough to 
alter the purely exponential heat decay. In other words, Eqs. (|2*7jl and ffi^i provide a demanding test of the coupling 
scheme under the present flow. Figure 9 shows the main Fourier component of the heat perturbation Q p obtained 
for 1- and 2-NHTcp. The dashed lines correspond to the theoretical expectation. It is evident that the l-NHT^p 
case does not obey the second law of thermodynamics at all. On the contrary, a rather good agreement with the 
theoretical trend is obtained when using at least 2-NHTcp. In Fig. 9, it is seen that the behaviour measured in the 

2- NHTcp case exhibits fluctuations around the theoretical straight line. Typically, the largest excursions last around 
3r from pure exponential decay. As previously stated, they may be due to the weakness of the entropy perturbation, 
but to confirm this statement we plan in the future to study some "heat-driven" flows (such as a heat pulse with 
initial zero mean velocity). 

VII. CONCLUSIONS 

We have presented the core of a hybrid continuum-particle method for moderate-to-large fluids which couples mass, 
momentum and energy transfers between two regions, C and P, described respectively by continuum fluid dynamics 
and by discrete particle Newtonian dynamics. Both domains overlap within a coupling region divided into two sub- 
cells which account for the two-way exchange: C^P and P^C. While the procedure at the P— >C cell is simply to 
average the particle-based (mass, momentum and energy) fluxes in order to supply open boundary conditions to the 
C domain, the operations at the C— >P cell are much less straightforward as they need to reconstruct a large number 
of (particles') degrees of freedom only from the knowledge of the three fluxes of conserved-quantities arising within C. 
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The present work has been concerned with extending the C— >P coupling to arbitrary rates of mass, momentum and 
energy transfer. To this end, the proposed method have been tested under unsteady flows which demand conformance 
to the whole set of conserved variable densities. In particular, we have considered the set of relaxing flows arising 
from hydrodynamics, namely longitudinal and transversal waves. We have followed the idea proposed by Flekkoy et 
al. in the sense that the scheme is explicitly based on direct flux exchange between the C and P regions. In order 
to deal with unsteady scenarios, we have shown that the fluxes injected into the particle region from the continuum 
region need to be measured exactly at the C— >P interface and not at the nodes of the continuum lattice. 

The implementation of flux exchanges requires the supply of energy currents to the particle system arising from 
the C domain due to advection, dissipation and conduction. To inject the correct amount of advected energy, the 
particle-averaged specific energy at the C^P cell needs to be equal to the continuum value. This can only be achieved 
if the new inserted particles are placed at positions where the (inter-particle) potential energy equals the C-specificd 
internal energy per unit mass. This severe condition has been implemented by the usher algorithm, whose purpose 
is two-fold: to provide the correct mass transfer rate and to ensure the balance of energy advection. In the proposed 
scheme, the balance of energy dissipation arises naturally, provided that the cell-averaged velocity and the injected 
momentum flux equal their continuum counterparts. This is made possible by applying the external force according 
to a flat distribution, instead of a biased one as used by Flekkoy et al. 5] but as a result, the new particles have 
to be inserted within a non- vanishing density environment. This is sorted out by the USHER in a very efficient way. 
Energy conduction has been implemented by using a set of Nose-Hoover thermostats adjacent to the C— +P interface, 
whose temperature and position are determined through the continuum local temperature gradient. Confirmation of 
the validity of this procedure is obtained from the correct rate of entropy production computed in our simulations of 
longitudinal waves. We showed that using only one thermostat per C^P cell (i.e. providing only the local value of T 
but not the heat flux) leads to negative entropy production. Therefore, in the context of energy transfer, this result 
reinforces the central importance of coupling-through- fluxes proposed by Flekkoy et al [j| . 

Enhancements to the present scheme are under investigation and merit some discussion here. The number density 
at C— >P may be controlled by a feedback algorithm which preserves the momentum flux balance. Other kinds of 
implementations for the energy transport by conduction also deserve to be considered. Finally, we plan to implement 
the P^C coupling in conjunction with a finite volume CFD solver in 3D. In order to extend the coupling scheme to 
higher dimensions some additional complications will need to be adressed. First, a mass flux will be assigned to each 
cell within an array of neighbourings C— >P cells. To adhere to mass continuity, particles will then have to be inserted 
within precisely defined finite regions and the insertion algorithm may have to pay an extra computational cost for 
this restriction of the search domain. We have checked, however, that the distance travelled by the usher algorithm 
from the initial trial position to the final insertion site is rather small (typically less than la and less than 0.5a on 
average) so we do not expect any significant extra cost if the search for insertion sites is done within volumes larger 
than (2a) 3 . In higher dimensions one may also have to smooth to some extent the variations of the mean mechanical 
quantities imposed along the C— >P region. To this end, it may be necessary to interpolate the external force along 
neighbouring C-^P cells in such a way that the local momentum flux imposed at each C— >P cell is still preserved. 
In the same way, although the Maxwell distribution, used here to choose the velocities of the new particles, proved 
to be sufficient to ensure momentum continuity for ID coupling (i.e, with no neighbouring C— >P cells), it may be 
convenient to use a Chapman-Enskog distribution in higher dimensions. This latter distribution enables the average 
velocity of the inserted particles to conform to the velocity gradient along neighbouring C— >P cells. We hope to report 
our findings in these areas in future publications. 
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overlapping region 




FIG. 1: (a) Spatial decomposition in our hybrid scheme. In this example the P region is adjacent to a physical surface 
represented by the rightmost shaded area. The continuum region spans the space to the left, at some distance from the surface. 
The overlapping region consists of a C— >P cell, where the C-flow is communicated to P, and a P^C cell, where particle 
averaged-fiuxes are injected into the C-flow. Dashed lines delimit the control cells of the C solver, with area A and grid-spacing 
AX. The letters O, W and E denote the centre of a cell and its west and east surfaces, respectively. The main cell's vectors 
(nw, n_E and npc) have been indicated (see text). 
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FIG. 2: Conditions imposed on the time step of the continuum solver Ate plotted versus the number density p. Variables are 
expressed in the LJ reduced units (a for length and r = (o- 2 m/e) 1//2 for time). As discussed in Sec. II V 21 Ate has to be greater 
than the collision time t co ; (the thick solid line) and smaller than AX/(2U) (indicated with dashed and dash-dotted lines). 
The typical flow velocity U is chosen to be either the sound velocity U = c s or the diffusive velocity U = v/L, according to 
each case discussed in Sec. II V 2l The grid-spacing is AX = 7r/(4fc ma x), where fc max is the largest wavenumber to be captured 
within the flow. 
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FIG. 3: The set-up for which the scheme has been tested (here region P is surrounded by C). The fluxes of continuum variables 
are imposed along the x direction, while y and z are periodic. The hybrid coupling is applied at the C— >P cells and the 
heat current is established along the Nose-Hoover thermostatted regions whose thickness is set between (3 — 4)o\ We used 
AX = (1 - 2)cr, L x = (10 - 40)cr and L y = L z = (7 - 9)<r. 



18 



(a) 



1 



:: •••««8#ooo oo0 ooo 00 

exp(-0.14t) 



A °0000o00 00o_ 



ooooo 




10 



FIG. 4: (a) The main Fourier component of the cell-averaged velocity v^l in (in units of ct/t, with r = (a 2 m/e) 1/ ^ 2 ). Results 
correspond to a transversal wave with wavenumber k — 0.31. Comparison is made between a calculation in which the full 
expression of the external force F ex was imposed (filled circles) and another which did not include the viscous contribution 
(open circles). The dashed line is the correct hydrodynamic decay, (b) The nondimensional autocorrelation function (ACF) of 
«^g in (i) for another transversal perturbation with k — 0.35, showing the theoretical decay in the (partially hidden) dashed line. 
In all cases the initial amplitude was v y ^l in = 1.0 and L x = 20a, L y = L z = la. In abscissas, time is nondimensionalized with 
the LJ reduced time unit r = (a 2 m/e) 1 ^ 2 . 
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FIG. 5: The total number of inserted particles J* 5N(t')dt' at the rightmost C— >P cell, along a simulation of a longitudinal wave 
with k — 0.168 inside a region (L x — 40a, L y = L z — 9a) with p e = 0.53 and T e = 3.5. The initial perturbative velocity profile 
was u x — 0.60 cos(fea:). The dashed line is the continuum prescription J* s(t')dt' (see text). Time has been nondimensionalized 
with r = (cr 2 m/e) 1/2 . 
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FIG. 6: The velocity at both C^P cells for the same parameters as in Fig. 5. The dashed line is the continuum prescription 
(partially hidden). In (a) we show the outcome of a calculation in which the momentum flux from the C-region was evaluated 
at x = xw (see Fig. 1); the instantaneous v x velocities are shown in lighter dotted lines, while thicker solid lines are the 
time-averaged velocities (v x ). In (b) we show the time-averaged velocities for another calculation with the same parameters as 
in (a), but evaluating the momentum flux at x — xo- All variables are nondimensionalised with the LJ potential units (it for 
length and r = ((j 2 m/e) 1 ^ 2 for time). 
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FIG. 7: (a) Time evolution of the Fourier transforms of a;-velocity i4 for the pure cosinusoidal perturbation of Fig. 5. In (b), the 
nondimensional autocorrelation function (ACF) of v£\i), corresponding to the pure cosinusoidal perturbation (solid line) and 
to another initial perturbation with |w^cos(0), v^ in (0) | = {0.60, 0.25} (dash-dotted line). In (a) and (b) the dashed line is the 

theoretical hydrodynamic solution. The remaining parameters are the same as those in Fig. 5. Variables are nondimensionalised 
with the LJ potential units (a for length and r = (cr 2 m/e) 1 / 2 for time). 
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FIG. 8: Time dependence of the Fourier amplitudes of the of (a) temperature T^, (b) density p^(= T, t exp(ikxi)/N), (c) 
energy per particle , and (d) pressure P^ ■ All quantities are nondimensionalised with the LJ potential units (a for length 
and e for energy). Comparison is made between a simulation with two Nose-Hoover thermostats per C^P cell (2-NHTcp), 
and another using only one (1-NHTcp, dotted line). In each figure, the analytical hydrodynamic expressions for the dominant 
Fourier component (in dashed lines) are explicitly written. The initial amplitudes were v^l OB (0) = 0.60, 2^(0) = —0.06, 
^(0) = 0.25, and p^(0) = 0.022. 
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FIG. 9: The main Fourier amplitude of the heat density perturbation — 



Q (1) 



(t) time-averaged along Ate = 1-0 (and 



multiplied by —1). The parameters are the same as those of Fig. 5 (and 8) and heat density is in units of e/a . The dashed 
line corresponds to the theoretical decay. Comparison is made between two (2-NHTcp) and one (1-NHTcp) Nose-Hoover 
thermostats per C^P cell. The latter violates the second law of thermodynamics. 



